//============================================================================
// Name        : csiena
// Author      : Josh Lospinoso
// Copyright   : MIT License
// Description : Austere, limited implementation of the R package RSiena in
//				 C++. Used for fast maximum likelihood and method of moments
//				 estimation of stochastic actor oriented models (SAOMs)
//				 without any frills or user experience niceties.
//============================================================================

#include "link.h"

namespace csiena {

ChainLink::~ChainLink() {
	delete(this->score);
	delete(this->information);
	delete(this->changeStats);
	delete(this->notionalScores);
}

ChainLink::ChainLink(int actors) {
	this->i=0;
	this->j=0;
	this->logProbability=0;
	this->score = new double[BETA_LENGTH];
	this->information  = new double[BETA_LENGTH*BETA_LENGTH];
	this->changeStats = new double[BETA_LENGTH*actors];
	this->notionalScores = new double[BETA_LENGTH*actors];
	for (int k=0; k<BETA_LENGTH; k++) {
		this->score[k] = 0;
		for (int l=0; l<BETA_LENGTH; l++) {
			this->information[k*BETA_LENGTH + l] = 0;
		}
		for (int j=0; j<actors; j++) {
			this->changeStats[BETA_LENGTH*j + k] = 0.0;
			this->notionalScores[BETA_LENGTH*j + k] = 0.0;
		}
	}

	this->actors = actors;
	this->previous=NULL;
	this->next=NULL;
}

ChainLink::ChainLink(const ChainLink& donorLink) {
	this->i=donorLink.i;
	this->j=donorLink.j;
	this->actors = donorLink.actors;
	this->logProbability=donorLink.logProbability;
	this->score = new double[BETA_LENGTH];
	this->information  = new double[BETA_LENGTH*BETA_LENGTH];
	this->changeStats = new double[BETA_LENGTH*this->actors];
	this->notionalScores = new double[BETA_LENGTH*this->actors];
	for (int k=0; k<BETA_LENGTH; k++) {
		this->score[k] = donorLink.score[k];
		for (int l=0; l<BETA_LENGTH; l++) {
			this->information[k*BETA_LENGTH + l] = donorLink.information[k*BETA_LENGTH + l];
		}
		for (int j=0; j<this->actors; j++) {
			this->changeStats[BETA_LENGTH*j + k] = donorLink.changeStats[BETA_LENGTH*j + k];
			this->notionalScores[BETA_LENGTH*j + k] = donorLink.notionalScores[BETA_LENGTH*j + k];
		}
	}

	this->actors = donorLink.actors;

	// Dangerous to copy these over, since links should never share the same neighbors.
	this->previous=NULL;
	this->next=NULL;
}

bool ChainLink::isLink() {
	return true;
}

double ChainLink::getLogProbability() {
	return this->logProbability;
}

double ChainLink::getScore(int betaIndex) {
	if (betaIndex < BETA_LENGTH) {
		return this->score[betaIndex];
	} else {
		std::cerr << "ChainLink::getScore: index out of bounds " << betaIndex << std::endl;
		return 0.0;
	}
}

double ChainLink::getInformation(int betaIndex1, int betaIndex2) {
	if ((betaIndex1 < BETA_LENGTH) & (betaIndex2 < BETA_LENGTH) ) {
		return this->information[betaIndex1 * BETA_LENGTH + betaIndex2];
	} else {
		std::cerr << "ChainLink::getInformation: index out of bounds " << betaIndex1 << "," << betaIndex2 << std::endl;
		return 0.0;
	}
}

double ChainLink::getChangeStatsForAlter(int betaIndex, int actorIndex) {
	if ((betaIndex < BETA_LENGTH) & (actorIndex < this->actors) ) {
		return this->changeStats[BETA_LENGTH*actorIndex + betaIndex];
	} else {
		std::cerr << "ChainLink::getChangeStatsForAlter: index out of bounds " << betaIndex << "," << actorIndex << std::endl;
		return 0.0;
	}
}

double ChainLink::getChangeStats(int k) {
	return this->getChangeStatsForAlter(k, this->j);
}


void ChainLink::evaluate(Digraph *x, Evaluation* fIn) {
	if (!fIn->isLinkEvaluation()) {
		std::cerr << "ChainLink::Evaluate error -- passed a CrimpEvaluation!" << std::endl;
	}
	LinkEvaluation* f = (LinkEvaluation*)fIn;
	f->preprocess(this->i, x);
	this->logProbability = log( f->getPij(this->j));

	// Calculate the "expected statistic", or the (weighted) average change statistic:
	double expectedStatistic[BETA_LENGTH];
	for (int k=0; k<BETA_LENGTH; k++) {
		expectedStatistic[k] = 0;
		for (int j=0; j< this->actors; j++ ) {
			expectedStatistic[k] += f->getPij(j) * f->getChangeStats(k,j) ;
			this->changeStats[BETA_LENGTH*j + k] = f->getChangeStats(k,j);
		}
	}
	// Calculate "notional scores," i.e. ignore which alter was actually selected,
	// and consider the score which would have resulted if the alter was chosen.
	// Also, take the opportunity to store the true score
	for (int k=0; k<BETA_LENGTH; k++) {
		for (int j=0; j< this->actors; j++ ) {
			this->notionalScores[k * BETA_LENGTH + j] = f->getChangeStats(k,j) - expectedStatistic[k];
		}
		this->score[k] = this->notionalScores[k * BETA_LENGTH + this->j];
	}
	// Calculate the information matrix from these ingredients
	for (int a=0; a<BETA_LENGTH; a++) {
		for (int b=0; b<BETA_LENGTH; b++) {
			this->information[a * BETA_LENGTH + b] = 0;
			for (int j=0; j<this->actors; j++) {
				this->information[a * BETA_LENGTH + b] += this->notionalScores[a * BETA_LENGTH + j] * f->getPij(j) * this->notionalScores[b * BETA_LENGTH + j];
			}
		}
	}
}

double ChainLink::evaluateProposal(Digraph *x, Evaluation* fIn) {
	if (!fIn->isLinkEvaluation()) {
		std::cerr << "ChainLink::evaluateProposal error -- passed a CrimpEvaluation!" << std::endl;
	}
	LinkEvaluation* f = (LinkEvaluation*)fIn;
	f->preprocess(this->i, x);
	return log( f->getPij(this->j));
}

void ChainLink::preprocess(Digraph *x, Evaluation* fIn) {
	if (!fIn->isLinkEvaluation()) {
		std::cerr << "ChainLink::evaluateProposal error -- passed a CrimpEvaluation!" << std::endl;
	}
	LinkEvaluation* f = (LinkEvaluation*)fIn;
	f->preprocess(this->i, x);
}

void ChainLink::print(int index) {
	printf(" (%3u) %3u %3u %1.3f ", index, this->i, this->j, this->logProbability);
	printf(" |");
	for (int k=0; k<BETA_LENGTH; k++) {
		printf(" %+3.2f", this->score[k]);
	}
	printf(" |");
	for (int k=0; k<BETA_LENGTH; k++) {
		printf(" %+3.2f", this->information[k * (BETA_LENGTH+1)]);
	}
	printf(" |");
	for (int k=0; k<BETA_LENGTH; k++) {
		printf(" %+3.2f", this->getChangeStats(k) );
	}
	printf("\n");
}

void ChainLink::simulate(Digraph *x, LinkEvaluation* f, MersenneTwister* r) {
	this->i = this->actors * r->next();
	f->preprocess(this->i, x);
	double temp = r->next();
	this->j = 0;

	while (temp > 0) {
		temp -= f->getPij(this->j);
		this->j++;
	}
	this->j--;
	this->evaluate(x, f);
}
}
